Multiscale or Multiresolution Optimization: An Application to 3D Gravity and Magnetic Data Inversion using Simulated Annealing

Presented by: Himanshu Bhardwaj, Founder & CEO, AICrux Business Solutions

Multiscale or multiresolution optimization is a powerful technique in numerical optimization used to tackle complex, high-dimensional problems. The idea is to progressively solve the problem at different levels of detail, starting with a coarse approximation and refining the solution at progressively finer scales. This approach is often employed in problems like image processing, machine learning, computational physics, and geophysical inversion where the dimensionality and complexity of the problem can make optimization challenging.

Key Theoretical Concepts

1. Dimensionality and Complexity:
   In high-dimensional optimization problems, finding the global optimum is computationally expensive and can be prone to getting stuck in local minima. By starting at a coarser resolution, the solution space is simplified, reducing the complexity of the optimization. Once a good approximation is found at this coarse scale, the solution is used as the starting point for higher resolutions.

2. Gradient Flow and Smoothness:
   At coarser scales, the objective function tends to be smoother, allowing faster convergence to a good approximate solution. This smoothness helps the optimization algorithm avoid small local minima, focusing on the global trends of the solution space. As the model is refined, the objective function becomes more detailed, and local minima may emerge, but the optimizer starts with a well-formed initial solution, making convergence faster.

3. Reduction of Search Space:
   The coarse model effectively reduces the search space, allowing the optimizer to explore fewer variables. As the resolution increases, the optimizer needs to work with more variables, but it does so with a better approximation of the solution from the previous coarse steps.

4. Hierarchical Model Refinement:
   In multiscale optimization, the model is refined step by step. Starting from a coarse approximation, the optimizer progressively works on finer representations of the same problem. The hierarchical nature allows the optimizer to balance global exploration at coarser scales and local exploitation at finer scales.

---

Steps for Implementing Multiscale Optimization

1. Coarse Initial Model (5x5x5 grid):

Start by creating a coarse approximation of the model, such as a 5x5x5 grid. This can be done by averaging or aggregating the values of multiple cells into a single value. For example, in a 3D grid, you could average groups of 2x2x2 cells from a finer resolution model to generate a coarser model. This step reduces the dimensionality and complexity, which speeds up the optimization process by focusing on larger-scale features of the solution.

Benefits: Reduces the number of variables and smooths the objective function, allowing faster convergence.

2. Convergence at the Coarse Level:

Apply the optimization algorithm (e.g., Simulated Annealing, Gradient Descent, or Genetic Algorithms) on the coarse model. This algorithm should run until the optimization reaches a point where the solution is considered satisfactory. The solution at this stage may not be the most precise but will provide a good approximation.

Key Point: The smoother objective function at the coarse scale helps the optimizer avoid local minima and focus on the global structure of the solution space.

3. Refining the Model:

Once the coarse model converges, refine it by increasing the resolution (e.g., from a 5x5x5 grid to a 10x10x10 grid). This step involves interpolating the coarser solution to create an initial guess for the next level of resolution. Interpolation techniques like linear interpolation, cubic splines, or bilinear interpolation are often used to project the coarse solution onto the finer grid.

Key Point: This refinement uses the solution from the coarser scale to initialize the finer-scale model, providing a good starting point for further optimization.

4. Optimize the Medium Model (10x10x10 grid):

With the refined model, rerun the optimization algorithm at the new resolution. The optimizer will now work on a more detailed version of the problem but with the advantage of having a good starting point from the previous coarse solution. This step allows the optimizer to explore finer details while still being guided by the global structure discovered at the coarser scale.

5. Final Fine Model (20x20x20 grid):

Finally, repeat the refinement process to reach the highest resolution (e.g., 20x20x20 grid). The optimization algorithm is then applied to this fine model, where it focuses on local details and further improves the solution. This final optimization step may take longer since the model is more complex, but the initial solution provided by the coarser models reduces the computational burden and helps guide the optimizer to the best possible solution.

Key Point: The process culminates in fine-tuning the solution at the highest resolution, where the optimizer balances global and local solution improvements.

---

Advantages of Multiscale Optimization

1. Efficiency in High Dimensions:
   One of the biggest challenges in high-dimensional optimization is the sheer size of the solution space. Multiscale optimization reduces this complexity by simplifying the problem at coarser levels. The algorithm converges faster and more efficiently than if it were applied directly at the highest resolution.

2. Better Convergence:
   Starting with a coarse model prevents the optimizer from getting stuck in local minima, which are more prevalent in high-resolution models. By progressively refining the model, the optimizer can explore both global and local minima more effectively.

3. Scalability:
   Multiscale methods are scalable to very large models, as the computational cost at each scale is relatively low. This makes them particularly useful for applications like computational fluid dynamics, geophysical inversion, and machine learning, where the problem size can be enormous.

4. Adaptability to Various Optimization Methods:
   This strategy can be adapted to many optimization algorithms, such as Simulated Annealing, Gradient Descent, Genetic Algorithms, and even Bayesian Optimization. It is flexible and can be tailored to specific problem domains by choosing appropriate coarsening and refinement techniques.

---

Practical Applications of Multiscale Optimization

1. Geophysical Inversion:

In geophysical studies, subsurface properties (such as density, velocity, or electrical conductivity) are estimated by inverting data from field measurements. Multiscale optimization helps refine models of the Earth's subsurface by starting with coarse approximations and progressively refining the resolution, leading to more accurate interpretations.

For example: In 3D Magnetic Data Inversion, accurately deciphering the 3D susceptibility model is crucial for identifying subsurface geological structures. By using multiscale optimization, the inversion starts with a coarse approximation of the subsurface susceptibility, allowing the algorithm to explore the broader features of the anomaly without getting trapped in local minima. As the optimization progresses through finer scales, the solution becomes more accurate, enabling better recovery of the spatial distribution of susceptibility. This multiresolution approach significantly enhances the clarity of subsurface features, making it particularly useful in magnetic data inversion, where high-dimensional models must be solved efficiently.

2. Image Processing and Computer Vision:

Multiscale techniques are widely used in image processing tasks like image segmentation and object detection, where coarse-to-fine strategies allow algorithms to focus on large, structural features before honing in on finer details.

3. Machine Learning and Neural Networks:

In deep learning, a coarse-to-fine strategy can be employed to improve the convergence of large neural networks. Coarse approximations of model parameters or reduced-resolution data can be used to speed up training before finer details are optimized.

4. Fluid Dynamics and Climate Modeling:

Simulations of fluid flow or climate systems often involve solving high-dimensional partial differential equations. Multiscale methods allow these simulations to start at coarser resolutions, making them computationally feasible while still capturing essential dynamics.

---

Conclusion

Multiscale or multiresolution optimization is an effective approach for solving high-dimensional optimization problems. By starting with a simplified, coarser version of the problem and progressively refining the model, it balances global exploration and local exploitation. This method improves convergence speed, reduces computational complexity, and provides a more robust solution by guiding the optimizer through a hierarchical refinement process.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import Forward_Model as fm
import time
from tqdm import tqdm
import grav_mag_inv as gmi

Let us create a 3D Forward Model synthetic data for Demonstration

# Parameters
L = 1000  # Size of the 3D grid (meters)
N = 20   # Number of points in each dimension
cylinder_radius = 200  # Radius of the cylinder (meters)
cylinder_length = 500  # Length of the cylinder along the Z-axis (meters)
cylinder_center = np.array([L/2, L/2])  # Center of the cylinder in the X-Y plane
susceptibility_cylinder = 0.04  # SI
susceptibility_background = 0.001  # SI
susceptibility, X, Y, Z = fm.create_3d_cylinder(cylinder_radius, cylinder_length, susceptibility_cylinder, susceptibility_background, L = L, N = N)

fig = fm.plot_3d_subsurface(susceptibility, X, Y, Z)
fig.show()

# Earth's magnetic field parameters
I = np.radians(49)  # Magnetic inclination in radians
D = np.radians(-10.2)  # Magnetic declination in radians
B_0 = 31000  # Earth's magnetic field strength in nT (assumed)

start_time = time.time()
delta_T_sim = fm.create_3d_forward_model(I, D, B_0, susceptibility, L, N)
elapsed_time = time.time() - start_time
print(f"Time taken: {elapsed_time:.2f} seconds")

Time taken: 0.78 seconds

fm.plot_forward_model(delta_T_sim, X, Y)

Inversion with coarse grid (4 X 4 X 4)

sus_coarse = gmi.aggregate_to_coarse(susceptibility, 5)

# Create a 3D mesh grid
x_5 = np.linspace(0, L, int(N/5))
y_5 = np.linspace(0, L, int(N/5))
z_5 = np.linspace(0, L, int(N/5))
X_5, Y_5, Z_5 = np.meshgrid(x_5, y_5, z_5, indexing='ij')

Down Sample the data

t_sim = gmi.downsample_observed_data(delta_T_sim, 5)

fm.plot_forward_model(t_sim, X_5, Y_5)

#X_init = susceptibility
X_init = susceptibility + np.random.uniform(0.0001, 0.001, size=8000).reshape(20,20,20)

sus_5 = gmi.aggregate_to_coarse(X_init, 5)

d_obs = t_sim
X_init = sus_5
N_5 = int(N/5)

best_X, best_error = gmi.simulated_annealing(I, D, B_0, d_obs, X_init, max_iters=10000, L = L, N = N_5, initial_temp=1.0, cooling_rate=0.95, lower_bound=0.001, upper_bound=0.9)
print("Best parameters found (20x20x20):", best_X)
print("Best error achieved:", best_error)

Iterations: 100%|███████████████████████████████████████████████████████████████| 10000/10000 [00:22<00:00, 439.71it/s]

Best parameters found (20x20x20): [[[0.00155824 0.00153119 0.00156794 0.0015479 ]
  [0.00151733 0.00155135 0.00155996 0.00156359]
  [0.00155983 0.00151801 0.0015305  0.00154988]
  [0.00151587 0.00156328 0.00154948 0.00151663]]

 [[0.00155749 0.00155419 0.00156715 0.00154228]
  [0.00153846 0.01869655 0.01868948 0.00158583]
  [0.00154447 0.01869018 0.0187067  0.00153029]
  [0.00154509 0.00157071 0.00154195 0.00157374]]

 [[0.00160322 0.00155371 0.00156697 0.00154237]
  [0.00155086 0.01876546 0.01871293 0.00152622]
  [0.00154637 0.01871566 0.01871384 0.00153085]
  [0.00156479 0.00152258 0.00155236 0.00156227]]

 [[0.00159152 0.0015635  0.00156933 0.00156427]
  [0.00155501 0.00158689 0.00155718 0.00154623]
  [0.00155106 0.00155195 0.00157276 0.00155689]
  [0.00160886 0.00156304 0.00154578 0.00157521]]]
Best error achieved: 378.3355228910471

Recovered Model

fig = gmi.plot_3d_subsurface(best_X, X_5, Y_5, Z_5)
fig.show()

start_time = time.time()
delta_T_sim_5 = fm.create_3d_forward_model(I, D, B_0, best_X, L, N_5)
elapsed_time = time.time() - start_time
print(f"Time taken: {elapsed_time:.2f} seconds")

Time taken: 0.00 seconds

Calculated Data

fm.plot_forward_model(delta_T_sim_5, X_5, Y_5)

Data Misfit

fm.plot_forward_model(t_sim - delta_T_sim_5, X_5, Y_5)

Inversion with Finer grid (8 X 8 X 8)

sus_20 = gmi.interpolate_solution(best_X, [20, 20, 20])

T_sim_8 = gmi.interpolate_solution(best_X, [8, 8, 8])
N_8 = int(N/2.5)

# Create a 3D mesh grid
x_8 = np.linspace(0, L, N_8)
y_8 = np.linspace(0, L, N_8)
z_8 = np.linspace(0, L, N_8)
X_8, Y_8, Z_8 = np.meshgrid(x_8, y_8, z_8, indexing='ij')

Start with the previous best parameters

fig = gmi.plot_3d_subsurface(T_sim_8, X_8, Y_8, Z_8)
fig.show()

T_obs_8 = gmi.downsample_observed_data(delta_T_sim, 2.5)

Observed Data (Downsampled 8 X 8 X 8)

fm.plot_forward_model(T_obs_8, X_8, Y_8)

d_obs = T_obs_8
X_init = T_sim_8

best_X, best_error = gmi.simulated_annealing(I, D, B_0, d_obs, X_init, max_iters=10000, L = L, N = N_8, initial_temp=1.0, cooling_rate=0.95, lower_bound=0.001, upper_bound=0.09)
print("Best parameters found (20x20x20):", best_X)
print("Best error achieved:", best_error)

Iterations: 100%|████████████████████████████████████████████████████████████████| 10000/10000 [02:18<00:00, 72.14it/s]

Best parameters found (20x20x20): [[[0.00155824 0.00154713 0.00153222 0.00153845 0.00155925 0.00156803
   0.00155699 0.0015479 ]
  [0.00154198 0.00153949 0.00153866 0.00154824 0.00156268 0.00156741
   0.00155925 0.00155293]
  [0.00151944 0.00152862 0.00154634 0.00155879 0.00156305 0.00156314
   0.00156218 0.00156157]
  [0.00152608 0.00153086 0.00153971 0.0015451  0.00154746 0.00155231
   0.00156034 0.00156469]
  [0.00155172 0.00154224 0.00152592 0.00151968 0.00152634 0.00153961
   0.00155279 0.00155877]
  [0.00155724 0.00154592 0.00152639 0.00151872 0.00152474 0.00153457
   0.00154151 0.00154405]
  [0.00153302 0.00153723 0.00154479 0.00154852 0.00154624 0.00153907
   0.00153015 0.00152573]
  [0.00151587 0.00153073 0.00155685 0.00156839 0.00156074 0.00154304
   0.00152494 0.00151663]]

 [[0.00155489 0.00154754 0.00153854 0.00154537 0.00156191 0.001567
   0.00155535 0.00154662]
  [0.00154242 0.00195346 0.00274097 0.00326903 0.00327985 0.00276389
   0.00197284 0.00155702]
  [0.00152495 0.00273102 0.00503302 0.00655455 0.00655693 0.00504756
   0.00276686 0.00157247]
  [0.00152941 0.00325192 0.00654014 0.00871463 0.00871752 0.00655442
   0.0032842  0.00157133]
  [0.00154857 0.00326023 0.00653034 0.00869927 0.00870804 0.00654781
   0.00327224 0.00155521]
  [0.00155348 0.00274454 0.00502066 0.00653179 0.00653852 0.0050302
   0.00274092 0.00154057]
  [0.00153671 0.00195477 0.00275131 0.00327404 0.00326847 0.00274089
   0.00194937 0.0015356 ]
  [0.0015246  0.00153972 0.00156443 0.00157026 0.00155696 0.00154242
   0.00153633 0.00153531]]

 [[0.00155447 0.00155195 0.00155016 0.00155664 0.00156609 0.00156509
   0.00155197 0.00154352]
  [0.00154752 0.00274735 0.00503874 0.00655636 0.00656028 0.00504881
   0.00276022 0.00156101]
  [0.00153669 0.00502996 0.01169678 0.01610085 0.01609759 0.01170117
   0.00506164 0.00158492]
  [0.00153569 0.00654086 0.0160953  0.02241234 0.02241279 0.01610487
   0.00657037 0.00157657]
  [0.0015429  0.00654264 0.01609031 0.02241189 0.0224215  0.01610853
   0.00655276 0.00154594]
  [0.00154692 0.00503351 0.01169222 0.01610246 0.01610994 0.01170311
   0.00503042 0.00153389]
  [0.00154368 0.00274983 0.0050498  0.00656411 0.0065556  0.00503607
   0.00274886 0.00155221]
  [0.0015407  0.00155337 0.00157147 0.00156805 0.00154941 0.00154297
   0.00155651 0.00156678]]

 [[0.00156975 0.00156435 0.00155692 0.00155923 0.00156661 0.00156403
   0.0015492  0.00153992]
  [0.00156201 0.00328178 0.00656496 0.00873627 0.00873551 0.00656183
   0.00327538 0.00155392]
  [0.0015478  0.0065595  0.01612309 0.02243762 0.02242523 0.01610557
   0.00656628 0.00157171]
  [0.00153972 0.00872464 0.02243784 0.03149863 0.03148849 0.02242382
   0.00873132 0.00156088]
  [0.00153999 0.00872025 0.02242876 0.03149708 0.03149905 0.02243095
   0.0087172  0.00153342]
  [0.00154416 0.00654921 0.01610654 0.02243316 0.02243958 0.01611514
   0.00654408 0.0015295 ]
  [0.00154847 0.00327063 0.0065576  0.0087296  0.00872773 0.00655617
   0.00327569 0.00155774]
  [0.00155044 0.0015538  0.00155716 0.00155168 0.00154343 0.0015478
   0.00156605 0.00157721]]

 [[0.0015932  0.0015798  0.00155838 0.00155451 0.00156415 0.00156441
   0.00154923 0.00153911]
  [0.00158001 0.00329746 0.00657651 0.00874569 0.00874335 0.00656393
   0.0032679  0.00154127]
  [0.00155555 0.00657495 0.01615162 0.02247044 0.0224493  0.01611004
   0.00654771 0.00154191]
  [0.00154101 0.00873813 0.02247152 0.03153762 0.03151445 0.02242758
   0.00871369 0.00153365]
  [0.00154036 0.00872961 0.02245229 0.03152227 0.03151314 0.02243213
   0.00871053 0.00152421]
  [0.00154647 0.00655242 0.01611099 0.02243719 0.02244187 0.0161162
   0.00654513 0.0015308 ]
  [0.00155324 0.00326727 0.00654235 0.00871571 0.00872556 0.00656051
   0.0032759  0.0015539 ]
  [0.00155641 0.00154736 0.00153281 0.0015303  0.00154073 0.00155432
   0.00156412 0.00156781]]

 [[0.00160448 0.00158776 0.00156026 0.001553   0.00156312 0.00156613
   0.00155413 0.00154555]
  [0.00158894 0.0027846  0.00506762 0.00657829 0.00657542 0.00505167
   0.00274665 0.00153905]
  [0.0015599  0.00506724 0.01175708 0.01616599 0.01614257 0.01170733
   0.00502483 0.00152762]
  [0.0015421  0.00657107 0.01616471 0.02249141 0.02246532 0.01611284
   0.00653556 0.00152261]
  [0.00154146 0.00656071 0.01614007 0.02246824 0.02245724 0.01611692
   0.00654152 0.00152732]
  [0.00155166 0.00504142 0.01170553 0.01611766 0.01612297 0.01171271
   0.00503737 0.0015398 ]
  [0.00156494 0.00275312 0.00502564 0.00653924 0.00655269 0.00504807
   0.00275704 0.0015548 ]
  [0.00157156 0.00155565 0.00153001 0.00152514 0.00154063 0.00155632
   0.00156178 0.00156215]]

 [[0.00159804 0.00158547 0.00156448 0.00155807 0.00156496 0.00156842
   0.00156236 0.00155761]
  [0.00158563 0.00199498 0.00277632 0.0032925  0.003289   0.00276238
   0.00196699 0.00155044]
  [0.0015614  0.00277387 0.00508371 0.00659877 0.00658007 0.00504458
   0.00274199 0.00153842]
  [0.00154394 0.00328154 0.00659406 0.00877279 0.00875585 0.00656115
   0.00326104 0.00153488]
  [0.00154185 0.00327118 0.00657221 0.00875437 0.00875293 0.00656954
   0.00326987 0.00154168]
  [0.0015571  0.00275401 0.0050416  0.00656119 0.0065698  0.00505654
   0.00275848 0.00155337]
  [0.00158234 0.00198334 0.00275105 0.00326426 0.00327019 0.00275806
   0.00197546 0.0015644 ]
  [0.00159585 0.00158118 0.00155444 0.00153986 0.00154197 0.00155261
   0.00156399 0.00156928]]

 [[0.00159152 0.00158237 0.00156693 0.00156175 0.00156639 0.00156958
   0.00156679 0.00156427]
  [0.00158168 0.0015795  0.00157482 0.00157022 0.00156639 0.00156284
   0.00155974 0.0015584 ]
  [0.00156156 0.00157077 0.00158428 0.00158268 0.00156753 0.00155377
   0.00154888 0.00154845]
  [0.00154496 0.00155673 0.0015762  0.00158137 0.00157105 0.00155737
   0.00154835 0.0015452 ]
  [0.00154187 0.00154672 0.0015573  0.0015677  0.00157294 0.00156873
   0.00155722 0.0015505 ]
  [0.00155964 0.00155547 0.00155061 0.00155528 0.00156594 0.00157028
   0.00156492 0.0015605 ]
  [0.00159149 0.0015812  0.00156215 0.00155094 0.00155118 0.00155812
   0.00156651 0.00157059]
  [0.00160886 0.00159638 0.00157113 0.0015506  0.00154286 0.00154975
   0.00156597 0.00157521]]]
Best error achieved: 290.5478258060534

fig = gmi.plot_3d_subsurface(best_X, X_8, Y_8, Z_8)
fig.show()

start_time = time.time()
delta_T_sim_8 = fm.create_3d_forward_model(I, D, B_0, best_X, L, N_8)
elapsed_time = time.time() - start_time
print(f"Time taken: {elapsed_time:.2f} seconds")

Time taken: 0.02 seconds

Calculated Data

fm.plot_forward_model(delta_T_sim_8, X_8, Y_8)

Data Misfit (with Finer grid 8 X 8 X 8)

fm.plot_forward_model(T_obs_8 - delta_T_sim_8, X_8, Y_8)

Inversion with finer grid (16 X 16 X 16)

Susc_16 = gmi.interpolate_solution(best_X, [16, 16, 16])
N_16 = int(N/1.25)

T_obs_16 = gmi.downsample_observed_data(delta_T_sim, 1.25)

d_obs = T_obs_16
X_init = Susc_16

# Create a 3D mesh grid
x_16 = np.linspace(0, L, N_16)
y_16 = np.linspace(0, L, N_16)
z_16 = np.linspace(0, L, N_16)
X_16, Y_16, Z_16 = np.meshgrid(x_16, y_16, z_16, indexing='ij')

best_X, best_error = gmi.simulated_annealing(I, D, B_0, d_obs, X_init, max_iters=1000, L = L, N = N_16, initial_temp=1.0, cooling_rate=0.95, lower_bound=0.001, upper_bound=0.09)
print("Best parameters found (20x20x20):", best_X)
print("Best error achieved:", best_error)

Iterations: 100%|██████████████████████████████████████████████████████████████████| 1000/1000 [02:42<00:00,  6.17it/s]

Best parameters found (20x20x20): [[[0.00155824 0.00155532 0.00154831 ... 0.00155604 0.00155032 0.0015479 ]
  [0.00155399 0.00155165 0.00154607 ... 0.00155671 0.00155141 0.00154917]
  [0.00154372 0.00154277 0.00154064 ... 0.00155831 0.00155413 0.00155237]
  ...
  [0.00153118 0.00153253 0.00153594 ... 0.00152898 0.00152587 0.00152468]
  [0.00152034 0.00152339 0.00153109 ... 0.00152534 0.00152064 0.00151886]
  [0.00151587 0.00151962 0.00152907 ... 0.00152396 0.00151863 0.00151663]]

 [[0.00155726 0.00155459 0.00154821 ... 0.00155564 0.00154998 0.0015476 ]
  [0.00155326 0.00155744 0.0015686  ... 0.00157905 0.00155771 0.00154925]
  [0.00154361 0.00156529 0.0016213  ... 0.00163909 0.00157744 0.00155333]
  ...
  [0.00153223 0.00155621 0.00161777 ... 0.00161135 0.00155091 0.00152743]
  [0.00152219 0.00153162 0.00155566 ... 0.00155061 0.0015307  0.00152302]
  [0.00151805 0.00152184 0.00153138 ... 0.00152704 0.0015229  0.00152137]]

 [[0.00155519 0.00155314 0.00154826 ... 0.0015546  0.00154908 0.00154677]
  [0.00155183 0.00157273 0.0016268  ... 0.00163644 0.00157382 0.00154931]
  [0.0015437  0.00162351 0.00182892 ... 0.00184654 0.00163711 0.00155554]
  ...
  [0.00153493 0.00161703 0.00182794 ... 0.00182299 0.00161513 0.00153431]
  [0.00152693 0.00155268 0.00161861 ... 0.00161545 0.00155634 0.00153345]
  [0.00152362 0.00152747 0.00153707 ... 0.0015348  0.00153359 0.00153323]]

 ...

 [[0.00159737 0.00159427 0.00158652 ... 0.00156239 0.00155955 0.00155836]
  [0.00159438 0.00161441 0.00166605 ... 0.00164085 0.00158021 0.00155656]
  [0.00158663 0.00166607 0.00187026 ... 0.00184271 0.00163348 0.00155207]
  ...
  [0.00158495 0.00166244 0.00186174 ... 0.00185099 0.00164561 0.00156564]
  [0.00159384 0.00161293 0.00166209 ... 0.00164505 0.00159015 0.00156875]
  [0.00159732 0.00159371 0.00158449 ... 0.00156486 0.00156851 0.00156991]]

 [[0.00159329 0.00159074 0.00158435 ... 0.00156536 0.00156342 0.0015626 ]
  [0.00159071 0.00159497 0.00160603 ... 0.00158604 0.00156805 0.00156101]
  [0.00158398 0.0016057  0.00166153 ... 0.00163938 0.00158008 0.00155702]
  ...
  [0.00159107 0.00161086 0.00166181 ... 0.00164653 0.00159112 0.00156955]
  [0.00160152 0.00160473 0.00161294 ... 0.00158887 0.0015771  0.00157252]
  [0.00160562 0.00160239 0.00159402 ... 0.00156637 0.00157158 0.00157364]]

 [[0.00159152 0.0015892  0.00158338 ... 0.00156653 0.00156495 0.00156427]
  [0.00158911 0.00158725 0.00158256 ... 0.00156474 0.00156336 0.00156278]
  [0.00158281 0.00158209 0.00158021 ... 0.00156034 0.0015594  0.00155905]
  ...
  [0.00159344 0.00159081 0.00158408 ... 0.00156701 0.00156999 0.00157113]
  [0.00160452 0.00160157 0.00159391 ... 0.00156706 0.00157211 0.0015741 ]
  [0.00160886 0.0016058  0.00159782 ... 0.001567   0.00157288 0.00157521]]]
Best error achieved: 400.25626973414313

fig = gmi.plot_3d_subsurface(best_X, X_16, Y_16, Z_16)
fig.show()

X_20 = gmi.interpolate_solution(best_X, [20, 20, 20])

start_time = time.time()
delta_T_sim_16 = fm.create_3d_forward_model(I, D, B_0, best_X, L, N_16)
elapsed_time = time.time() - start_time
print(f"Time taken: {elapsed_time:.2f} seconds")

Time taken: 0.16 seconds

Calculated or Predicted Magnetic Data (16x16x16)

fm.plot_forward_model(delta_T_sim_16, X_16, Y_16)

start_time = time.time()
delta_T_sim_20 = fm.create_3d_forward_model(I, D, B_0, X_20, L, N)
elapsed_time = time.time() - start_time
print(f"Time taken: {elapsed_time:.2f} seconds")

Time taken: 0.46 seconds

fm.plot_forward_model(delta_T_sim_20, X, Y)

import GravMagPro as gm

True Model XZ slice at Y = 500 m

col = gm.get_color()
plt.contourf(X[:,10,:], Z[:,10,:], susceptibility[:,10,:], cmap=col)
plt.colorbar()
plt.title('Original Susceptibility with Cylindrical Anomaly')
plt.xlabel('X (m)')
plt.ylabel('Z (m)')
plt.ylim(Z[:,10,:].max(), Z[:,10,:].min())
plt.show()

Recovered Model XZ slice at Y = 500 m

plt.contourf(X[:,10,:], Z[:,10,:], X_20[:,10,:], cmap=col)
plt.colorbar()
plt.title('Recovered Susceptibility with Cylindrical Anomaly')
plt.xlabel('X (m)')
plt.ylabel('Z (m)')
plt.ylim(Z[:,10,:].max(), Z[:,10,:].min())
plt.show()

3D Plots of True and Recovered Model

fig = fm.plot_3d_subsurface(susceptibility, X, Y, Z)
fig.show()

fig = gmi.plot_3d_subsurface(X_20, X, Y, Z)
fig.show()